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ABSTRACT 

Aims. We investigate the stability of non-isothermal Bonnor-Ebert spheres in the context of a model that includes a self-consistent 
calculation of the gas temperature. In this way, we can discard the assumption of equality between the dust and gas temperatures, and 
study the stability as the gas temperature changes with the chemical evolution of the cooling species. 

Methods. We use a gas-grain chemical model to calculate the chemical evolution. The model includes a time-dependent treatment 
of depletion onto grain surfaces, which strongly influences the gas temperature as the main coolant molecule CO depletes from the 
gas. The dust and gas temperatures are solved with radiative transfer calculations. For consistent comparison with previous work, we 
assume that the cores are deeply embedded in a larger external structure, corresponding to visual extinction Ay Xt = 10 mag at the core 
edge. We also study the effect of lower values of Ay 1 . 

Results. We find that the critical non-dimensional radius fi, determining the maximal density contrast between the core center and 
the outer boundary, derived here is similar to our previous work where we assumed Tdust - 7g as ; the gi values lie below the isothermal 
critical value ~ 6.45, but the difference is less than 10%. We find that chemical evolution does not affect notably the stability 
condition of low-mass cores (< 0.75 M e ) which have high average densities and a strong gas-grain thermal coupling. In contrast, for 
higher masses the decrease of cooling owing to CO depletion causes substantial temporal changes in the temperature and in the density 
profiles of the cores. In the mass range 1-2 M s , decreases with chemical evolution, whereas above 3 M 0 , instead increases with 
chemical evolution. We also find that decreasing A“* strongly increases the gas temperature especially when the gas is chemically old, 
and this causes £i to increase with respect to models with higher Ay 1 . However, the derived £i values are still close to £o- The density 
contrast between the core center and edge derived here varies between 8 to 16 depending on core mass and the chemical age of the 
gas, compared to the constant value ~ 14.1 for the isothermal BES. 
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1. Introduction 


The Bonnor-Ebert sphere ( Bonno r 195 6. lEbert 1955 c hereafter 
BES), i.e., an isothermal gas sphere in hydrostatic equilibrium, 
has been used succesf ully to approximate the density struc¬ 
tures of prestellar cores ( Bacmanne^alJ l200(il : I Alves et al 1200 it 
iKandori et all 120051: iMarsh et al.1 2014) . H owever, the assump¬ 
tion of isothermicity is not va l id generally (IZucconi et al .11200 It 


Ward^homjrson^eyil. 200$, iPagani et al.1 20041 Craosi et all 


20071: Juvela & Ysardl 201 ll) . To accommodate for a radial tem¬ 


perature profile, a non-isothermal version of the BES (referred 
here to as a modified Bo nnor-Ebert sphere; MBES) has been 
studied in the lite rature (Evans et al.l l200lb iGalli et all 120021 
iKeto & Fieldll2005l: ISinila et al.ll201 lh. 


In previous studies of the MBES, it is either assumed that the 
dust and gas temperatures are equal, or the gas temperature has 
been derived based on some standard abundances for the cool¬ 
ing species. However, both approaches are approximations and 
do not hold generally. In Sipila et al. (2011 : hereafter Paper I), 
we studied the stability of MBESs that are deeply embedded in 
a larger external structure, e.g., a molecular cloud, correspond¬ 
ing to a high visual extinction Ay 1 = 10 mag at the edges of the 
studied cores. We also assumed that T & ust = 7 gils , which holds 
well for low-mass cores with high average densities, but is not 
valid at lower density where the collisional coupling between 


gas and dust is weak. In the present paper, we aim to generalize 
the analysis of Paper I by including in the stability calculations 
a self-consistent determination of the gas temperature. This is 
accomplished by calculating chemical evolution in model cores 
with a com prehensive gas-grain chemical model (Sipilal 12012 : 
ISipila et all 1 201 3l) . followed by a determination of the gas tem¬ 
perature with a radiative transfer model at different time steps, 
taking advantage of the time-dependent chemical abundances. In 
this way, we can study the stability condition not only as a func¬ 
tion of the core mass but also as a function of “chemical time”, 
which is here defined as the chemical age of the core since some 
initial state (see Sect. 12.21) . We can determine whether the results 
of Paper I are significantly affected when Tj ust A T gas , especially 
in high-mass cases when the cores have low average densities 
and consequently weak gas-grain thermal coupling. 

The paper is organized as follows. In Sect. [2) we discuss our 
model calculations in detail. In Sect.[3] we present the results 
of our calculations, and discuss them in Sect. [4] We present our 
conclusions in Sect.0 


2. Method 

This section outlines how the calculations are carried out in prac¬ 
tice. The chemical model is discussed here only briefly, and we 
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refer the reader to ISipilal d2012l) and ISinila et aP (l2013h for a 
complete description of the model. 


2.1. The MBES 


In what follows, we discuss the basic properties of the MBES in 
a rather concise form; A more detailed discussion on the MBES 
can be found in Paper I. 

The MBES differs from the BES in that it is non-isothermal. 
Assuming the ideal gas equation of state and hydrostatic equi¬ 
librium, and that T - T(r), the density distribution of the MBES 
is given by 

1 d I r 2 dp dT 
r 1 dr\p dr ^ dr 

This equation can be transformed into non-dimensional form by 
making the substitutions 


) 4nGmp 


P = 




( 2 ) 


r=/3 l/2 A~ l/2 {. 


(3) 


In the above, A is the central density of the core; t(£) = T (£)/'!£, 
where T c is the central temperature of the core; /3 = kTJ4nGm, 
where k and m are the Boltzmann constant and the average 
molecular mass of the gas (assumed here equal to 2.33 amu), re¬ 
spectively; £ and il/(£) are dimensionless variables (£ represents a 
non-dimensional radius). With these substitutions, Eq. iQ} trans¬ 
forms to the modified Lane-Emden equation 


,_2 d_ pi ^ 
d£ d£ 



(4) 


We impose the following boundary conditionsat the core center: 
ijj — 0, difr/d£ — 0, r — 1 and dr/di' = 0. These boundary 
conditions ensure that T — T c and p — A at £ = 0, i.e., at the 
center of the core. To solve Eq. 0, one also needs to supply 
a temperature profile, which has to be determined externally - 
we discuss the temperature calculations in more detail below. 
Finally, the density profile of the MBES is given by substituting 
the solution function i/» and the temperature profile into Eq. 0. 


Fix core mass M and non-dimensional radius £ oul 




Construct a BES corresponding to the chosen values of M and 6u, 


1 


Iteration 



Calculate a dust temperature profile 


1 


Solve Eq. 0 


1 


Calculate central density from Eq. 0 


t 


Calculate density profile from Eq. 0 


1 


Divide the core into concentric shells and calculate 
chemistry separately in each shell 


Iteration 



t 


Calculate the gas temperature profile 


I 

Solve Eq. 0 

\ 


Calculate central density from Eq. 0 


1 


Calculate density profile from Eq. 0 


1 


^ Calculate a new dust temperature profile 

± 

Extract fi, A and t//(£ 0 u t ) 


Repeat the above for different values of 


Examine the stability of the solutions using Eq. 0 


2.2. Determining the stability of the MBES 

The non-dimensional radius £ is a free parameter - for a given 
core mass, there exist a series of core configurations correspond¬ 
ing to different values of £. In what follows, the non-dimensional 
radius of each core configuration will be represented by £ out . It 
should be noted that the solution to Eq. 0) is unique to each core 
configuration (see Paper I), so that the solution is applicable in 
the interval 0 < € < fout for each value of fou. (i.e., separately 
for each core configuration). 

As discussed in Paper I, determining the stability of an 
MBES of given mass is analogous to finding the core config¬ 
uration which is critically stable. We construct the different core 
configurations in the same way as in Paper I, i.e., by following an 
iterative process. To illustrate the process, outlined in Fig.[T] let 
us fix the core mass and the non-dimensional radius f out . We first 
construct a BES corresponding to the chosen values of core mass 
and ( ; out , and adopting T — 10 K. After this step, we determine 
a dust te mperature profile for the BES using r adiative transfer 
modeling (IJuvela & Padoani [20031; ljuvelali2005b . The dust tem¬ 
perature profile is then used to solve Eq. 0. I n Pape r I, we con¬ 
sidered two different dust models (from lOssenkonf & Henningl 


Fig. 1. Steps taken in deriving the critically stable configuration 
of an MBES of given mass. See text for detailed explanations of 
each step. 

fl994l and [Li & Draineii2001l) . but here we only consider the for¬ 
mer. This issue is discussed in Sect. 14. II 

The solution to Eq. 0 and the dust temperature profile are 
used to derive the central density of the MBES from the equation 

r«(^) (5) 

where r out = T(£ oat )/T c . After this step, the density profile of 
the MBES is calculated from Eq. 0. To ensure the consistency 
of the density profile with the temperature profile, we perform 
additional iterative calculations, calculating a new dust tempera¬ 
ture profile and a new density profile in sequence. Both profiles 
typically reach their final solutions in a few iterations. 

Next, we proceed to calculate the gas temperature. In Paper I, 
we made the simplifying assumption that T & ust = 7 gas . However, 
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this is a rather crude estimate for cores with total hydrogen den¬ 
sity nu ^ 10 5 cm 3 , where the coupling between the dust and the 
gas is weak. In this paper, we carry out an updated analysis of 
the results of Paper I, replacing the dust temperature with the 
gas temperature in all of the relevant formulae. The ga s tem¬ 
peratu re is solved identically to the method discussed in ISipilal 
( 2012! ). That is, we divide the first-approximation MBES den¬ 
sity profile into concentric shells, solve the chemical evolution 
separately in each shell and extract chemical abundance profiles 
for the cooling species (see below) as functions of radial dis¬ 
tance from theatre center. The chemical model is adopted from 
ISinila et all (l2013b ; the model includes a full description of gas- 
grain chemistry and the deuterated forms of chemical species 
with up to 4 atoms. 

The gas temperature i s calc ulated with a Monte Carlo ra¬ 
diative transfer program dJuvelal 1 997 ) which balances heat¬ 
ing and cooling functions to solve the gas temperature. In the 
present model, the gas is heated externally by cosmic rays and 
by the photoelectric effect (lGoldsmithll2001 ). However, the lat¬ 
ter process is only important at Ay < 2 mag (Liuvela & Ysardl 
1201 ll) . Molecular line cooling is calculated for the following 
species: 12 CO, 13 CO, C ls O, C, O and O 2 . Since the adopted 
chemical model does not include the isotopes of the various 
species, we adopt the isotopic ratios 12 CO/ 13 CO = 60 and 
12 CO/C ls O = 500 (IWilson & Roodl 11994 . We also include 
the energy exchange (collisional coupling) between the gas a nd 
the dust grains following the description of lGoldsmithl (1200 llPl 
We implicitly assume that chemical timescales are longer than 
the radiative and dynamical timescales, so that a core is able to 
quickly reach a new equilibrium state as the chemistry evolves. 

After the gas temperature profile has been calculated, it is 
used to solve Eq. i(4j again, producing a new density profile 
which is then used to redetermine successively the dust tem¬ 
perature, the chemical abundances and the gas temperature as 
outlined above. The iteration is repeated a few times until the 
density and the (gas and dust) temperature profiles converge to¬ 
ward their respective final solutions. 

Carrying out the above iterative process for different values 
of^oun while keeping the core mass constant, we obtain p, A and 
i/r as functions of foot and use these to construct the pressure 
derivative 


2 p [ 

p-'dp + /I” 1 (5/1 - | 

(?) 

1 , &4out 

out 

3V 

[p- l * 5p- A~ l 6A + 2^ 

rat &Tout] 


see Paper I for details on how this expression is derived. Finally, 
the critically stable configuration corresponds to the lowest value 
of £out for which the pressure derivative is zero - all values_of C„n 
abo ve this value correspond to unstable configurations dBonnorl 
1 956L Paper I). In the next section, we present the critical &ut 
values for a range of core masses derived as outlined above, and 
we label the critical £ out value for each core mass as 4 : i . Similarly, 
we label the values of p, A and |A(fout) corresponding to critical 
cores as P\, A\ and respectively. 

Chemical abundances vary with time, and this influences 
the gas temperature as 1) atomic species are processed into 
molecules and 2) the main coolant species CO in its various 
isotopic forms freezes onto the dust grains. To investigate the 
effect of varying chemical abundances on the stability of the 

1 We note that I Young et all d2004f) have considered a stronger gas- 

grain coupling than in the model of Goldsmith d200lh . which might 

somewhat compensate for the depletion of coolant molecules at inter¬ 

mediate densities if adopted in our model. 


MBES, we have carried out the above analysis for four differ¬ 
ent time steps, corresponding to 5 x 10 4 , 1 x 10 5 , 5 x 10 5 or 
1 x 10 6 years of chemical evolution since the initial state of the 
core (see below). In practice, we calculate £1 for each core mass 
four times, extracting the chemical abundances at either 5 x 10 4 , 
1 x 10 5 , 5 x 10 5 or 1 x 10 6 years. The different gas temperature 
profiles at the four time steps translate to marked differences in 
density profiles, and also in the critical 4 ; i values. This issue is 
discussed in the next section. We note that we assume external 
visual extinction of A“ l = 10 mag in the chemical and radiative 
transfer calculations, so that the model cores are assumed to ex¬ 
ist deeply embedded inside larger structures, such as molecular 
clouds. This choice ensures the possibility of direct comparison 
of our results with those of Paper I. 

In the chemical calculations, we assume that the gas is ini¬ 
tially atomic, with the exception of hydrogen which is in molec¬ 
ular form (!Sipilall2012l) . The choice of the chemical composition 
of the gas in the beginning of the calculation (labeled as to in the 
figures) is a parameter of the model. The four timesteps defined 
above measure the extent of chemical development with respect 
to the atomic initial state. We note that an atomic chemical com¬ 
position may not be consistent with an initial physical structure 
corresponding to a BES (as is assumed here) because chemi¬ 
cal processing of the gas is expected to take place during the 
formation of the BES itself. However, here we do not attempt 
to impose constraints on “absolute” core ages, but to investi¬ 
gate if chemical evolution can influence the stability condition. 
Therefore, the chosen four timesteps are simply representative 
points along the chemical development track of a given core and 
its surroundings. 

3. Results 

We present in Fig. [2] the critical non-dimensional radius 4 ; i and 
other related quantities (see below) as functions of core mass and 
chemical time, as derived from Eq. (|6j. The values of 4’i derived 
here are consistently lower than in Paper I for the OH94 model 
(their Fig. 4). However, the 4'i values predicted by the two works 
agree to within 10% regardless of core mass. Hence we derive 
similar stability for (deeply embedded) critical MBESs regard¬ 
less of whether we assume 7d U st = T gas or 7d ust 4 r gas . 

Evidently, there is marked variation in 4u both as a function 
of core mass and as a function of chemical time. A clear de¬ 
creasing trend in £1 is seen from the lowest masses (~ 0.1 M 0 ) 
up to about 0.75 M 0 . 4 : i coincides with the isothermal critical 
radius at ~ 0.25 M B . For masses above 0.75 M 0 , 4u is, at early 
times, constant to an accuracy of ~ 0.15 regardless of core mass. 
However, at late times the value of £1 depends on the core mass. 
In the range 0.75 M Q - 2.5 M 0 , £1 decreases as a function of 
chemical time while for M > 2.5 M Q an increase with chemi¬ 
cal time is evident. For an isothermal BE sphere, the translation 
from 4”out to physical quantities is straightforward: a larger ^ Q ut 
translates directly to a higher density contrast between the cen¬ 
ter and the edge, while the physical radius depends on 4'om. the 
central density, and the isothermal temperature. Here, the con¬ 
version between the various parameters is less evident because 
the temperature is also a function of radius, and the temperature 
profile affects not only the density contrast (Eq.[2} but also the 
solution to the Lane-Emden equation and the determination of 
the central density (Eqs.[4]and[5}. To understand the significance 
of the temporal changes in £ 1 , we have to consider the associated 
changes of A\, P\ (i.e., T c j), and ij/\. 

The central density of a critically stable non-isothermal 
sphere, Ai, decreases monotonously with increasing mass like 


3 







































O. Sipila et al.: On the stability of non-isothermal Bonnor-Ebert spheres. II. The effect of gas temperature on the stability 








Fig. 2. Critical non-dimensional radius £) \ (uppei left~htmd panel), i/4 £) i) — ijj \ (upper middle panel), density contrast dj/yciout.i — 
(T’outa/T’c.i)^ 1 (upper right-hand panel), central temperature T c \ - 4nGm[3\/k (lower left), central density A\/m (lower middle), 
and the outer core radius R out ,\ = [i i:1 A (lower right) of each critical core configuration as functions of core mass. The profiles 
are plotted at four different time steps, given in the legend in the bottom middle panel. The blue dotted lines in the upper panels 
represent the corresponding critical values for the (isothermal) Bonnor-Ebert sphere (£ : o = 6.45, i//q = 2.64, and Ao/p ou t.o = 14.1). 


in the isothermal case. Chemical evolution starts to modify A\ 
above 0.75 M Q where it increases with chemical time. The cen¬ 
tral temperature T c first increases as a function of mass, reaches a 
maximum, and then decreases again towards the largest masses. 
Like Ai, T c shows hardly any temporal variation below 0.75 M Q . 
Beyond this point, however, T c increases strongly with chemical 
time and the maximum shifts towards larger masses. Finally, iff\ 
decreases with chemical time for masses below ~ 2.5 M Q , but 
increases for the higher masses. 

The approximate constancy of i;\, A\ and fi\ over chemical 
time at the lowest masses (< 0.75 M Q ) implies that the outer 
radii of critically stable low-mass cores are similar regardless 
of how long the gas has been chemically processed. For higher 
core masses, the physical size of a critical core always decreases 
with chemical time, but the effect is much more prominent for 
the highest masses. This is a consequence of the increasing aver¬ 
age temperature caused by chemical evolution; hotter cores are 
better able to withstand gravity and external pressure, and thus 
higher densities are required for the cores to collapse. 

For critically stable low-mass MBESs, the density con- 
strast is always lower than that of a critical BES (~ 14). For 
M < 0.75 M q , the density contrast increases with chemical 
time owing to changes in the temperature profile; the tempera¬ 
ture at the edge of the cores increases with CO depletion (see 
also Sect. 14.21 note that the central temperature is not affected 
owing to the efficient gas-grain thermal coupling at the core cen¬ 
ter) and the consequent increase in the T out /T c ratio overcompen¬ 
sates the decrease of e' ! ‘ due to the decrease of Aj as a function of 


chemical time. For M > 0.75 M Q , the density contrast largely 
follows the changes in i/y. For M > 2.5 Mq, the density contrast 
can slightly exceed that of the critical BES (~ 14). 

The described tendencies are controlled by three pro¬ 
cesses: 1) thermal coupling between gas and dust, 2) molecu¬ 
lar line cooling, and 3) depletion of molecules onto dust grains. 
The gas-grain thermal coupling is efficient at densities above 
/;(Hi) ~ 10 5 cm~ 3 . Therefore the lowest-mass stable cores 

which have the highest (central and average) densities are prac¬ 
tically unaffected by changes in the chemical composition of the 
gas, because the temperature profile is determined by the inter¬ 
action between the dust and the interstellar radiation field. The 
same is true for the centers of slightly more massive cores up to 
about 1 M q . However, their outer parts are cooled efficiently by 
molecular line emission, until molecules, in particular CO, start 
to freeze out. With the diminished cooling, the gas temperature 
rises. 

The most massive MBESs have low central densities and 
the temperature is determined by dust heating and molecular 
line cooling throughout. The temperature increases strongly in 
their central parts owing to CO depletion. For each time step, 
there is a local maximum in the central temperature (as a func¬ 
tion of core mass) which shifts toward higher core masses for 
longer chemical times. This is because higher-mass cores have 
increasingly lower average densities and hence the gas-grain 
thermal coupling is progressively weaker, allowing the temper¬ 
ature to increase further with CO depletion. However, the de¬ 
pletion timescales are very long in the outer layers where the 


4 













O. Sipila et al.: On the stability of non-isothermal Bonnor-Ebert spheres. II. The effect of gas temperature on the stability 



Fig. 3. Density (black lines, scale on left y-axis) and gas temperature (blue lines, scale on right y-axis) profiles of an MBES with 
M — 1 M 0 , £ out = 5 as functions of core radius. The profiles are plotted for four different time steps, labeled in the left panel. Left 
panel: External Ay = 10 mag. Middle panel: External Ay = 2 mag. Right panel: External Ay = 1 mag. 


densities are low, and line cooling can operate there also at late 
chemical times. The steep temperature gradient results in a large 
density contrast between the center and edge for the most mas¬ 
sive cores (~ 16 for the critical cores). 

Given the relatively straightforward temporal changes in T c 
and Ti, we can deduce that the function t//| is mainly responsible 
for the mass-dependent behavior of f i as a function of chemical 
time. The changes in the cooling efficiency alter the radial den¬ 
sity distribution described by i/y, which leads either to a decrease 
or an increase of depending on the core mass. 


4. Discussion 

4.1. The effect of temperature on the results 

We have presented the critical radii of MBESs that are deeply 
embedded in a parent molecular cloud, corresponding to a high 
external visual extinction Ay Xt = 10 mag. This choice facilitates 
comparison with the results of Paper I where the same assump¬ 
tion was made. However, as pointed out in Sect. [2 the determi¬ 
nation of the temperature may have a strong impact on our re¬ 
sults because the solution to Eq. i(4} is not unique. To study how 
our results might vary with different initial assumptions, we have 
calculated the properties of a 1.0 M e core with two low values 
of A“ 1 (2 mag and 1 mag). 

Figure [3] presents the density and temperature profiles of a 
1.0 M 0 core with (arbitrarily chosen) £out = 5 at different time 
steps, assuming Ay 1 = 10, 2 or 1 mag. The gas temperatures 
are similar at the core center regardless of the choice of ex¬ 
ternal Ay, because of the high density: even for A“ l = 1, the 
visual extinction at the core center is ~ 10 mag (we assume 
n(H 2 )/Ay = 9.4 x 10 2 °: lBohlin et al.Hl97& ). However, there are 
significant differences between the gas temperature profiles out¬ 
side the core center depending on the choice of external Ay. This 
is because a decrease of visual extinction increases the heating 
by the photoelectric effect in the core. Although the photodis¬ 
sociation rate of the major coolant molecule CO is strongly in¬ 
creased as well, the net effect on core cooling is not large be¬ 
cause in these circumstances atomic carbon takes over as the 
main coolant species. The gas temperature at the edge of the 
core for A™ 1 = 1 and t - 10 6 yr is ~ 26 K. 



Fig. 4. Critical non-dimensional radius (upper left panel) of 
the 1.0 Mq MBES as a function of chemical time for different 
values of Ay', indicated in the figure. 


Evidently, the central densities are similar in all cases, even 
though the determination of the central density depends also on 
the temperature and <// profiles, and not only on the central tem¬ 
perature (Eq.0. The density profiles are also similar for chemi¬ 
cal times up to t = to + 5 x 10 s yr. For t = to + 1 X 10 6 yr at low 
Ay 1 , the strongly increased gas temperature changes the mass 
distribution in the core so that the slope of the density profile 
becomes somewhat steeper. 

The differences in temperature depending on the choice of 
Ay* translate to very different solutions to Eq. ©, and conse¬ 
quently the non-dimensional radius of the critical configuration 
is different in all cases. Figure[4]shows the change of the critical 
non-dimensional radius £\ of the 1.0 M e core as a function of 
chemical time, assuming Ay 1 = 10, 2 or 1 mag. Evidently, the 
critical radius increases as A“' decreases. Also, for low Ay 1 , the 


5 
















O. Sipila et al.: On the stability of non-isothermal Bonnor-Ebert spheres. II. The effect of gas temperature on the stability 





R [AU] xlO 4 R [AU] xlO 4 R [AU] xlO 4 

Fig. 5. Thermal pressure ( left-hand panel), density (middle panel ) and gas temperature ( right-hand panel) of an M = 1.0 M 0 MBES 
with £,ut = 6 as functions of distance from the core center at four different time steps, labeled in the left-hand panel. The vertical 
lines in each panel represent the outer radius of the core at the different time steps. 


critical radius tends to increase as a function of chemical time, 
whereas for Ay 1 = 10 mag, there is a slight decreasing trend. 

We find that the f) values are within ~ 15 % of each other 
regardless of the choice of Ay 4 , when M - 1.0 M 0 . For lower 
core masses, we expect the results to be closer together because 
of the high average densities, which lead to significant temper¬ 
ature differences only at the very edge of the cores. However, 
for M > 1.0 M 0 , the average density is small and hence the Ay 
gradient through the core is shallow, and we expect temperature 
effects to be more prominent. The increase of with decreasing 
Ay implies that stable isolated cores can have a larger density 
constrast (between the core center and the edge) than cores em¬ 
bedded in molecular clouds. 

We note that in Paper I, we found similar values for re¬ 
gardless of the du st m odel, although the dust temperature given 
by the model of U_&Drain^j 2001k is low er than that given by 
the model of lOssenkonf & Hennindd 19941) (see Pape r I). In the 
present paper, we have not considered the [Li & Drainel d2001l) 
dust model, as we do not expect our results to be strongly de¬ 
pendent on the dust model. This is motivated on the one hand by 
the results of Paper I and on the other hand by the fact that 
does not change significantly even if the gas temperature profile 
is radically different (Fig.[4]i. 


4.2. Core evolution 

In the stability calculations presented in this paper, we consid¬ 
ered a series of model cores defined by the mass and the non- 
dimensional radius, so that we could derive values for the crit¬ 
ical non-dimensional radius according to an analytical formula 
(Eq.[6]t. However, the choice of the parameters used to define and 
MBES is free and when constructing a model core one could, in¬ 
stead of the non-dimensional radius, set for example the external 
pressure or the (dimensional) outer radius. In a realistic scenario, 
one can expect all of these quantities to change as the core and its 
surroundings evolve; regardless of how the two required param¬ 
eters are chosen, implicit assumptions are made simultaneously 
on the properties of the medium outside the core. 

If we consider a fixed non-dimensional radius like earlier in 
this paper, the properties of the core are limited by the external 
thermal pressure exerted on it. This can be understood by study¬ 
ing the evolution of the medium at the core boundary. In Fig.|3 
we plot the thermal pressure, density and gas temperature pro¬ 


files of an MBES with M - 1.0 M Q and £ out = 6. The £ grid 
has been converted to physical radius R using Eq. <|3}, and verti¬ 
cal lines marking ( ; out = 6 have been inserted at each timestep to 
help estimate the change in R. The gas temperature increases as a 
function of time owing to CO depletion. For late chemical times 
with increasing thermal pressures, the configuration defined by 
ftut = 6 decreases in size (and increases in avera ge de nsity) . 
We note that the hydrodynamical models of lKeto & Field! (120051 1 
predict that MBESs can exhibit oscillatory behavior. Expansion 
as a function of time is not achieved in our models because the 
gas temperature is nearly always increasing, leading to more ef¬ 
ficient compression of the cores at long chemical timescales. 

The changes in physical radius, density profile etc. at the dif¬ 
ferent time steps considered here do not necessarily represent 
evolutionary tracks for the cores, because in a realistic scenario 
neither the non-dimensional radius nor the mass of a core are 
constant. However, the present study demonstrates that the con¬ 
ditions for the core stability change with chemical evolution. For 
example, in the case of the 1 M 0 MBES, we expect a core con¬ 
figuration with fout ~ 6 to be critically stable if the medium is 
chemically young, while the same configuration would be un¬ 
stable if associated with chemically old gas. Nevertheless, when 
a core has been found, e.g., based on the column density and 
temperature distributions, to agree with a MBES model, the di¬ 
agrams presented in Fig.[2]of this paper can be used to estimate 
its stability. 


5. Conclusions 

We analyzed the stability of non-isothermal (modified) Bonnor- 
Ebert spheres with a new model that includes a self-consistent 
determination of the gas te mperat ure. We compared our results 
with those of ISinila et al.l < 201 lb . where it was assumed that 
Tdust = Tgas- We found that the critical non-dimensional radius 
changes with the chemical evolution in the core and its sur¬ 
roundings, especially for cores with M > 1 M Q . This is because 
of the depletion of the coolant species (mainly CO) onto grain 
surfaces which raises the gas temperature. Therefore, cores that 
exist in a chemically young environment are expected to have 
a different stability condition than those that exist in regions of 
chemically old gas. 

The i : i values derived in the present work are slightly lower 
than in the models of lSinila et al. ! ( 2011 ). However, the results of 
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the two works agree to within 10 %. The bulk of our analysis was 
carried out for deeply embedded cores with external visual ex¬ 
tinction A™ 1 = 10 mag; test calculations for lower values of Ay 1 
yield higher critical radii. In summary, our results indicate that 
the stability of the modified Bonnor-Ebert sphere is similar re¬ 
gardless of whether one assumes 7’ c i ust = r gas or 7’ ( | llst + 7' gas and, 
by extension, also similar to th e stabi lity of the classic isothermal 
Bonnor-Ebert sphere (see iSipila et al.l20TTh - when one assumes 
that the core is deeply embedded in a larger structure. 
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